Optimization of Johnson–Cook Constitutive Model Parameters Using the Nesterov Gradient-Descent Method

Numerical simulation of impact and shock-wave interactions of deformable solids is an urgent problem. The key to the adequacy and accuracy of simulation is the material model that links the yield strength with accumulated plastic strain, strain rate, and temperature. A material model often used in engineering applications is the empirical Johnson–Cook (JC) model. However, an increase in the impact velocity complicates the choice of the model constants to reach agreement between numerical and experimental data. This paper presents a method for the selection of the JC model constants using an optimization algorithm based on the Nesterov gradient-descent method. A solution quality function is proposed to estimate the deviation of calculations from experimental data and to determine the optimum JC model parameters. Numerical calculations of the Taylor rod-on-anvil impact test were performed for cylindrical copper specimens. The numerical simulation performed with the optimized JC model parameters was in good agreement with the experimental data received by the authors of this paper and with the literature data. The accuracy of simulation depends on the experimental data used. For all considered experiments, the calculation accuracy (solution quality) increased by 10%. This method, developed for selecting optimized material model constants, may be useful for other models, regardless of the numerical code used for high-velocity impact simulations.


Introduction
Numerical simulation of high-velocity interactions of deformable solids requires using material models that adequately describe material behavior at high strain-rates and temperatures. The role of shock-wave processes increases with an increase in the impact velocity. Known material models, such as the classical empirical Johnson-Cook (JC) model [1,2], that describe the change in yield strength (both hardening and softening) give different results compared to the experimental data with increasing plastic strain, plastic strain rate, and temperature. Classical material models developed for calculating the reliability of construction elements and structures at high impact velocities can give critically incorrect results. An adequate numerical description of the behavior of materials under dynamic loading leads to modifications in known empirical models, increasing the number of fitting parameters and therefore increasing the complexity of such models. Most researchers use commercial software products, in particular ANSYS/LS Dyna, that implement known material models with standard constants, causing difficulties with the adequate simulation of high-velocity shock-wave processes.
The JC constitutive model is widely used in many scientific and engineering applications. For example, Zhang C. et al. [3] used the JC model to simulate the deformation and fracture of fiber-reinforced polymer composites at normal and elevated temperatures. Xie H. et al. [4] calibrated the JC model parameters for melt-cast explosives under parameters to ensure good agreement of numerical calculations with experimental data. This trend has no tendency to decrease. This paper presents the numerical calculations of classical Taylor rod-on-anvil impact tests for cylindrical OFHC copper samples. A method for selecting the JC model constants using an optimization algorithm based on the Nesterov gradient-descent method [42] was proposed. To optimize the selection of the JC model constants, a solution quality function was proposed, which can estimate the deviation of the calculations from the experimental data and determine the optimal JC model parameters. The optimal JC model parameters were selected.

Formulation of the Problem
A 3D simulation of a cylindrical projectile impact with a non-deformable target (Taylor test) was performed. The material of the cylindrical sample (oxygen-free high-conductivity copper, OFHC) was chosen to compare the numerical results with the experimental data [41]. Simulations of the experiments [16,43] were conducted for copper specimens with corresponding initial conditions (geometric dimensions, velocity).
The 3D calculations were performed using an elastic-plastic medium model and our own research software (COMP3, v1.0) based on the modified finite element method developed by G.R. Johnson [44]. The system of basic equations and relations of the finite element method can be found in [45][46][47]. The COMP3 software has shown its adequacy for a wide range of high-velocity interactions of deformable bodies [48,49].
The method used 4-node tetrahedral finite elements. A finite-element mesh for a cylinder was generated as follows. The cylinder was divided along the height by planes parallel to the end surfaces into N Z -1 layers of the same height, where N Z is the number of nodes along the vertical Z axis. Each layer was divided into N Y -1 equal sectors, where N Y is the number of rays, and into N R -1 rings of equal thickness, where N R is the number of nodes along the ray. As a result, the cylinder was divided into quadrangular and triangular (near the axis) prisms. The quadrangular prisms consisted of 6 tetrahedral elements, while the triangular prisms consisted of 3. Figure 1 shows the finite element model of a cylinder. This problem is symmetric, so 1/2 of the cylinder is modelled. The number of elements in the finite element model was 1314, and the number of nodes was 6120 (N R = 9, N Y = 9, N Z = 18). parameters to ensure good agreement of numerical calculations with experimental data. This trend has no tendency to decrease. This paper presents the numerical calculations of classical Taylor rod-on-anvil impact tests for cylindrical OFHC copper samples. A method for selecting the JC model constants using an optimization algorithm based on the Nesterov gradient-descent method [42] was proposed. To optimize the selection of the JC model constants, a solution quality function was proposed, which can estimate the deviation of the calculations from the experimental data and determine the optimal JC model parameters. The optimal JC model parameters were selected.

Formulation of the Problem
A 3D simulation of a cylindrical projectile impact with a non-deformable target (Taylor test) was performed. The material of the cylindrical sample (oxygen-free high-conductivity copper, OFHC) was chosen to compare the numerical results with the experimental data [41]. Simulations of the experiments [16,43] were conducted for copper specimens with corresponding initial conditions (geometric dimensions, velocity).
The 3D calculations were performed using an elastic-plastic medium model and our own research software (COMP3, v1.0) based on the modified finite element method developed by G.R. Johnson [44]. The system of basic equations and relations of the finite element method can be found in [45][46][47]. The COMP3 software has shown its adequacy for a wide range of high-velocity interactions of deformable bodies [48,49].
The method used 4-node tetrahedral finite elements. A finite-element mesh for a cylinder was generated as follows. The cylinder was divided along the height by planes parallel to the end surfaces into NZ-1 layers of the same height, where NZ is the number of nodes along the vertical Z axis. Each layer was divided into NY-1 equal sectors, where NY is the number of rays, and into NR-1 rings of equal thickness, where NR is the number of nodes along the ray. As a result, the cylinder was divided into quadrangular and triangular (near the axis) prisms. The quadrangular prisms consisted of 6 tetrahedral elements, while the triangular prisms consisted of 3. Figure 1 shows the finite element model of a cylinder. This problem is symmetric, so 1/2 of the cylinder is modelled. The number of elements in the finite element model was 1314, and the number of nodes was 6120 (NR = 9, NY = 9, NZ = 18).

Modification of the JC Constitutive Model
The JC constitutive model [1,2] was used for the 3D simulation of deformation of a copper cylinder. The disadvantage of this model is that an increase in plastic deformation with an increase in the strain rate leads to an unlimited increase in the yield strength, which is not observed in experiments. In particular, Zelepugin S.A. et al. [41] showed that OFHC copper in high-velocity impact experiments with plastic strain accumulation demonstrates an increase in microhardness from 900 MPa to~3000 MPa. Moreover, the JC constitutive model parameters were experimentally determined by the model developers at relatively low strain-rates (up to 10 3 s −1 ). However, the strain rates reached values of 10 6 s −1 and more at high-velocity impacts. Model constants determined for numerical simulation of high-velocity processes at low strain-rates led to discrepancies between numerical and experimental results. This raised the question of determining the parameters of the JC constitutive model for high strain-rates.
Numerical calculations using the original JC model at an impact velocity of 316 m/s revealed a hardening of ≥5 times due to an increase in equivalent plastic strains, which is 2 times higher than the hardening measured in [41]. In this regard, the JC constitutive model was modified by limiting the maximum value of plastic hardening in terms of the dimensionless factor B max : Here, σ 0 is the quasi-static yield strength (denoted as A in the original model [1]); ε pl is the equivalent plastic strain; .
ε 0 is the dimensionless plastic strain rate; . ε pl is the plastic strain rate; . ε 0 is the initial strain rate ( . ε 0 = 1 s −1 ); T 0 , T, and T m are the initial, current, and melting temperatures, respectively; and B, n, C, and m are the material model constants.

Solution-Quality Function
We proposed to introduce a solution quality function to quantitatively estimate the discrepancy between experimental and numerical results.
The reliable measured parameters of the cylinder after impact are the residual length of the cylinder L and the radius of points on the cylinder lateral surface W f ( Figure 2). A less-reliable measured parameter is the maximum cylinder radius, R f , due to the possible fracture of the cylinder in contact with the rigid wall. It is worth noting that the specimens were partially broken along the external radius at impact velocities of 316 m/s and above [41].
For more accurate agreement between experimental data and numerical results, we proposed to estimate five points on the lateral surface of the cylinder. The solution quality function Q f was selected in the form as follows: where The weighting factors 20, 2, and 1 were selected in such a way that the residual length, being a reliable measured parameter, had the greatest effect on the solution quality. Deviations of the lateral surface radius had less effect on the solution quality. The smallest contribution was given by the maximum radius, as it varied in the experiments during the tension of the external edge of the cylinder. Note that the quality function could be chosen in different ways, and the influence of quality function was not investigated in this paper. For more accurate agreement between experimental data and numerical results, we proposed to estimate five points on the lateral surface of the cylinder. The solution quality function f Q was selected in the form as follows: The weighting factors 20, 2, and 1 were selected in such a way that the residual length, being a reliable measured parameter, had the greatest effect on the solution quality. Deviations of the lateral surface radius had less effect on the solution quality. The smallest contribution was given by the maximum radius, as it varied in the experiments during the tension of the external edge of the cylinder. Note that the quality function could be chosen in different ways, and the influence of quality function was not investigated in this paper. The Q f function can numerically estimate the deviation of calculations from the experimental data and find the model parameters using optimization methods. Q f is a positive value and turns to zero when the calculations of the external surface of the cylinder after impact match the experimental data. The external surface of the cylinder remains smooth after deformation, so the chosen number and location of radius measuring points W f is adequate to approximate the shape of the external surface with high accuracy.
The introduced quality function Q f solves the problem of finding model parameters that minimize this function and are optimal.
Since the number of parameters to be optimized was small (6: σ 0 , B, n, B max , C, m), gradient-descent methods were effective for solving this problem. However, it should be noted that calculating the Q f value for a single numerical experiment takes a long time (up to several hours, Intel Xeon 3.2 GHz). Therefore, the number of calculation steps should be minimized. For this purpose, the Nesterov gradient-descent method was chosen to determine the Q f minimum [42].
The idea of the Nesterov gradient-descent method is that the direction of the gradient is calculated at the predicted value of the currently found point, and the movement is performed relative to this point. The convergence rate of the method increases if the gradient of the function changes lightly. In this case, the distance between the points grows, and the speed of movement along the gradient increases. At the same time, sharp changes in the gradient direction with error lead to a reduction in the movement step, and the method effectively finds the local minimum.
Numerical calculations showed that the influence of parameters m and n on the numerical results was weak for the considered impact conditions. In addition, the parameter σ 0 is a quasi-static yield strength and cannot change at high strain-rates from a physical point of view. Thus, three parameters of the modified model remained under consideration: B, B max , and C. The JC model parameters were optimized to determine the minimum of the function:

Numerical Results and Discussion
The parameters of the Johnson-Cook model for OFHC copper were estimated using the experimental data of the Taylor test [41] and well-known experiments on the impact of cylinders with a rigid wall by Gust W.H. [16] and by Wilkins M.L. and Guinan M.W. [43]. An extended review of the calculated and experimental data for different medium models is given in [50].
For simulation, we chose experiments with the initial data shown in Table 1. The first series of calculations was performed using the original Johnson-Cook model parameters for copper, shown in Table 2. in the gradient direction with error lead to a reduction in the movement step, and the method effectively finds the local minimum. Numerical calculations showed that the influence of parameters m and n on the numerical results was weak for the considered impact conditions. In addition, the parameter σ0 is a quasi-static yield strength and cannot change at high strain-rates from a physical point of view. Thus, three parameters of the modified model remained under consideration: B, Bmax, and C.
The JC model parameters were optimized to determine the minimum of the function:

Numerical Results and Discussion
The parameters of the Johnson-Cook model for OFHC copper were estimated using the experimental data of the Taylor test [41] and well-known experiments on the impact of cylinders with a rigid wall by Gust W.H. [16] and by Wilkins M.L. and Guinan M.W. [43]. An extended review of the calculated and experimental data for different medium models is given in [50].
For simulation, we chose experiments with the initial data shown in Table 1. The first series of calculations was performed using the original Johnson-Cook model parameters for copper, shown in Table 2.   Figure 3a shows the fields and isolines of the specific shear strain energy, and Figure  3b shows the temperature fields and isolines at times of 30 and 87 µs. It can be seen that the maximum shear strains, as well as temperature, were generated in the center of the cylindrical specimen at the contact boundary between the specimen and the rigid wall. The maximum temperature reached 949 K by 87 µs. The calculations were terminated at 87 µs when the average velocity of the specimen along the Z axis became positive. We did not use a failure criterion and friction in our simulation [51].   Figure 3a shows the fields and isolines of the specific shear strain energy, and Figure 3b shows the temperature fields and isolines at times of 30 and 87 µs. It can be seen that the maximum shear strains, as well as temperature, were generated in the center of the cylindrical specimen at the contact boundary between the specimen and the rigid wall. The maximum temperature reached 949 K by 87 µs. The calculations were terminated at 87 µs when the average velocity of the specimen along the Z axis became positive. We did not use a failure criterion and friction in our simulation [51].

(m/s) T (K) Reference
Calculated and experimental profiles of the external surfaces of the cylinders are shown in Figure 4.
The values of the solution quality function are presented in Table 3.     Table 2); 1, 2, 3, 4, 5, 6 are the sequence number of tests in Table 1.
The values of the solution quality function are presented in Table 3. Calculations showed that the original JC model constants did not agree quantitatively with the experimental data. Tests 3 and 4 demonstrated a near-perfect match between numerical and experimental data for the cylinder's final length. However, the shape of the lateral surface and the maximum diameter differed. For test 1, the calculation described well the experimental shape of the lateral surface, but there was a discrepancy in the cylinder's final length. Other tests also showed a discrepancy between numerical and experimental results. It is worth noting that the simulation performed for tests 1 and 2 were in good agreement with the calculations presented in [50].
Optimization of the JC model parameters was performed as follows: 1. Optimization of parameters for each of tests 1, 2, 3, 4, 5, and 6; 2. Optimization of parameters for tests 1 and 2: ( ) Materials 2023, 16, x FOR PEER REVIEW 9 of 13 Figure 5 shows the calculated and experimental profiles of the external surface of the cylinders when optimizing the parameters for each of experiments 1, 2, 3, 4, 5, and 6. The JC model parameters after optimization are presented in Table 4.  Table 4) for tests 1-6. In tests 1-4, the hardening limit was reached in a small volume and did not affect the process. The material deformed heavily in tests 5 and 6, so the results were affected by the limitation of the maximum plastic hardening. For this reason, in tests 5 and 6, the numerical optimization reduced the hardening limit Bmax to a value of ~2.8, which corresponds to the microhardness experiments [41]. However, in the experiments [41], the increase in maximum microhardness at individual points reached 3.5 times, and the average microhardness increase at the front end of the cylinder increased by 1.8-2.3 times.  Table 4) for tests 1-6. In tests 1-4, the hardening limit was reached in a small volume and did not affect the process. The material deformed heavily in tests 5 and 6, so the results were affected by the limitation of the maximum plastic hardening. For this reason, in tests 5 and 6, the numerical optimization reduced the hardening limit B max to a value of~2.8, which corresponds to the microhardness experiments [41]. However, in the experiments [41], the increase in maximum microhardness at individual points reached 3.5 times, and the average microhardness increase at the front end of the cylinder increased by 1.8-2.3 times. Table 5 compares the calculated optimal model parameters with the original parameters [2] for different tests, where B/B 0 and C/C 0 are the ratios between optimal and original parameters, Q f0 is the solution quality before the optimization of parameters, and Q f is the solution quality after the optimization of parameters. Good agreement between the calculated and experimental data was observed at Q f < 0.07 (tests 1, 2, 5, 1 + 2, 3 + 4 + 5 + 6). It should be noted that the original parameters of the Johnson-Cooke model better described the behavior of copper in the experiments in [41] at impact velocities of 162 m/s and 167 m/s than in the experiments in [16,43]. However, at higher impact velocities, the deviation of the optimal model parameters from those determined in the low-velocity experiments increased up to two times.
Calculations showed that the experiments in [16,43] are best described by the Johnson-Cooke model, assuming that the OFHC copper parameters differ from the measured parameters [2] (B-decrease by 30%, C-decrease by 10%). At impact velocities of 162 m/s and 167 m/s [41], the experiments are described exactly by the same set of parameters. However, with an increase in the impact velocity, the deviation of the calculated values from the experimental values increases. The set of experimental data [41] is better described under the assumption that hardening parameter B of the used copper increases by a factor of 1.9, while parameter C decreases slightly.
Determining a set of parameters that optimally describes the set of experiments considered [16,41,43] showed that the average solution quality improves by 10%, but this improvement is achieved by the accurate description of the experiments in [16,43] (by 1.7 and 1.8 times) and a slight deviation (by 20% and 10%) in the description of the experiments in [41] for impact velocities of 225 m/s and 316 m/s, respectively.

1.
The JC constitutive model was modified by introducing a material-hardening limit for plastic deformation, B max , at high strain-rates.

2.
A solution quality function, Q f , was proposed to estimate the deviation of calculations from the experimental data. The final length of the cylinder, the radius of the lateral surface of the cylinder at five points, and the maximum radius of the cylinder were taken as the function parameters, with weighting factors of 20, 2, and 1 according to the effect on the final quality of the solution and reliability of the parameter measurement.

3.
An optimization algorithm for selecting parameters B and C of the JC constitutive model and the limiter B max was developed to find the best agreement between the calculated and experimental data for the Taylor impact test using the Nesterov gradientdescent method. 4.
The optimal parameters, namely, B, B max , and C, of the modified JC constitutive model were calculated for nine sets of experimental data. The solution quality in some experiments increased by several times when using optimal parameters. For all experiments, the solution quality improved by 10% after optimization. 5.
The developed method for optimizing the selection of the constitutive model constants can be adapted for a wide range of problems (arbitrary set of optimized parameters, arbitrary material models, and software codes, including ANSYS/LS Dyna). Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.